/*
 * GraphAlgorithm.h
 *
 *  Created on: 02/09/2012
 *      Author: marcio
 */

#ifndef GRAPHALGORITHM_H_
#define GRAPHALGORITHM_H_

#include "DirectlyAddressedSet.h"
#include "BitMapDirectlyAddressedSet.h"
#include "Graph.h"
#include "heap.h"
#include <cstring>

/*!
 * \class ColoringHeuristic
 *
 * \brief This class encapsulates a set of algorithms to color a graph
 *
 * The most basic ones are to color or to determine the connected
 * components of a graph.
 *
 * A coloring of a graph is a assignment of colors to the vertexes of the graph in such away
 * that if to vertexes are adjacent they have different colors.
 *
 * The connected components of a graph are a set of graphs such that each vertex and edge belong
 * to only one graph, the union of its vertexes sets and edge sets are the vertex set and the edge
 * set of the initial graph and all the graphs in the set are connected, what means that in each graph
 * there are a path between every pair of vertices.
 *
 * \tparam T define the bitmap directly addressed sets
 *
 * \see Graph
 * \see BitMapDirectlyAddressedSet
 *
 */
template <typename T, template<typename T> class D>
class ColoringHeuristic {
protected:

	const long			n;			//!< max number of elements

	Factory<D<T> > * 	factory;
	D<T> * const 		bitR;
	D<T> ** const 		colors;
	AndFilter<T> * 		andFilter;

	long *				invR;		//!< inverse of input mappings
	D<T> * const 		degR;

	/*!
	 * \brief Creates a covering algorithm object
	 *
	 * Free list is initialized containing all nodes.
	 *
	 * \param n size of internal arrays and sets
	 * \param factory a factory to BitMapDirectlyAddressedSet to be used to instantiate the variables.
	 *
	 */
	ColoringHeuristic(long n, Factory<D<T> > * f);

public:
	virtual ~ColoringHeuristic();

	/*!
	 *
	 *
	 * \param g the input graph
	 * \param res the vector in which the colors are returned: the position \c i in this vector indicates the color
	 * of the vertex \c R[i] in the greedy coloring of the subgraph induced by \c R.
	 * \param R a vector representing the verticees to be colored and the order to be used to color them.
	 * \param r the number of verticees in the subgraph.
	 *
	 * \return the number of colors used in the greedy coloring of the subgraph induced by \c R returned in \c cover
	 *
	 */
	long run(Graph<D<T> > * g, long * res, long * R, long r) = 0;
	static long bucketSort(long * color, long cmin, long cmax, long * R, long r);

	/*!
	 * \brief Made a greedy coloring of the subgraph induced by the vertexes in R, with the order that they
	 * appear in R.
	 *
	 * Made a coloring, greedy coloring of the subgraph induced by the vertexes in R. The order used to
	 * built suc coloring is the same order that the vertexes appear in R. The algorithm visit each element
	 * in R, in order, and look for the smallest color not used in one of its neighbours to color it. If
	 * there is not such color a new color will be add to the set of possible colors.
	 *
	 * \param n size of internal arrays and sets
	 * \param factory a factory to BitMapDirectlyAddressedSet to be used to instantiate the variables.
	 *
	 * \return the covering algorithm object created
	 */
	static ColoringHeuristic<T, D> * greedy(long n, Factory<D> * factory);

	/*!
	 * \brief Determines a greedy coloring of the subgraph induced by the specified vertices using the Dsatur algorithm.
	 *
	 * The order used to guide the greedy coloring is specified by a function given as input.
	 * The Dsatur algorithm keeps track of the so called
	 * saturation degree of each vertex. The saturation degree of a vertex \c v is the number of already colored vertices
	 * in the neighborhood of \c v.
	 * At each iteration, the uncolored vertex with the largest saturation degree (with the largest degree used to break ties, whenever
	 * possible) is chosen to be colored.
	 *
	 * \param n size of internal arrays and sets
	 * \param factory a factory to BitMapDirectlyAddressedSet to be used to instantiate the variables.
	 *
	 * \return the covering algorithm object created
	 *
	 * \see http://www.math.tu-clausthal.de/Arbeitsgruppen/Diskrete-Optimierung/publications/2002/gca.pdf
	 */
	static ColoringHeuristic<T, D> * dsatur(long n, Factory<D> * factory);

	/*!
	 * \brief Cover the vertexes of the graph according its degree and the degree of its neighbours.
	 *
	 * Return an order to the vertexes in R, according the degrees in the subgraph induced for it.
	 * More clearly, if cover[i] = j it means that in the graph induced by the vertexes R[k] that has
	 * cover[k] < j, R[i] is a vertex with minimum degree.
	 *
	 * The integer value returned by this function is the biggest k such that k is the minimum degree
	 * in the subgraph induced by R[0] .... R[i] for all i < r.
	 *
	 * \param g the graph
	 *
	 * \param cover represent the order over the vertexes of g, the position i of this vector
	 * indicate the position of the vertex i in the order.
	 *
	 * \param R a vector representing the vertexes to be colored and the order to be used to color then.
	 *
	 * \param r the number of vertexes in the subgraph.
	 *
	 * \param factory a factory to OrderedSet to be used to instance the variables.
	 *
	 * \return the maximum mininum degree in the subgraph.
	 */
	static ColoringHeuristic<T, D> * mcr(long n, Factory<D> * factory);

	/*!
	 * \brief
	 *
	 * \param g
	 *
	 * \param k
	 *
	 * \param rank
	 *
	 * \return
	 *
	 */
	static double FCP2(Graph * g, long k, long** rank);
};

template<typename T, template<typename T> class D>
ColoringHeuristic<T, D>::ColoringHeuristic(long n, Factory<D<T> > * f) : n(n), factory(f) {
	bitR = factory->newInstance(n);
	colors = new T[n];
	long j;
	for (j = 0 ; j < n; j++)
		colors[j] = factory->newInstance(n);
	andFilter = new AndFilter<T>();
	degR = f->newInstance(n);
	invR = new long[n];
}

template<typename T, template<typename T> class D>
ColoringHeuristic<T, D>::~ColoringHeuristic() {
	delete andFilter;

	for (long j = 0 ; j < n; j++)
		delete colors[j];
	delete[] colors;
	delete bitR;
	delete degR;
	delete[] invR;
}

template <typename T, template<typename T> class D>
class TripleList {
protected:

	long *				next1;		//!< list of keys, level 1
	long *				next2;		//!< list of keys, level 2
	long *				prev1;		//!< list of keys, level 1
	long *				prev2;		//!< list of keys, level 2
	long *				key1;		//!< key, level 1
	long *				vert;		//!< first element of the associated list of level 3
	long *				node1;		//!< node in the list of level 1

	long				freehead;	//!< free nodes of the lists of degrees
	long				tail12;		//!< lists of levels 1 and 2
	long				head1;		//!< head of the list of level 1

	long *				next3;		//!< list of keys, level 3
	long *				prev3;		//!< list of keys, level 3
	long * 				key3;		//!< key, level 3
	long * 				node2;		//!< node in the list of level 2

	long				tail3;		//!< lists of level 3

	TripleList(long n, Factory<D<T> > * f);
	~TripleList();

	long removeFreeHead();
	void removeFree(long j);
	void addFreeHead(long i);
	long insertLevel1Before(long k);
	long removeLevel1(long j);
	long removeLevel2(long j);
	long removeLevel3(long i);

	/*!
	 * \brief Inserts the specified vertex in the list of structural degrees
	 *
	 * The insertion is performed in a node *after* the specified basepoint.
	 * This method searches for the node in the list of structural degrees until a member corresponding to the specified node's degree
	 * is found or it is detected that such a member does not exists.
	 *
	 * \param j the basepoint
	 * \param i the vertex to insert
	 *
	 * \return the node of the list of structural degrees that receives the specified vertex
	 */
	long insertLevel2(long j, long i);
};

template<typename T, template<typename T> class D>
TripleList<T, D>::TripleList(long n, Factory<D<T> > * f) {
	long nn = 1 + (n << 1);  // an index is left as sentinel

	next1 = new long[(nn << 2) + (nn << 1) + 1];
	next2 = next1+nn;
	prev1 = next2+nn;
	key1 = prev1+nn;
	next3 = key1+nn;
	prev3 = next3+n+1;
	key3 = prev3+n+1;
	node2 = key3+n+1;

	prev2 = prev1;
	node1 = key1;
	vert = next1;

	freehead = 0;
	for (tail12 = 0; tail12 < nn; tail12++) {
		next1[tail12] = tail12+1;
		prev1[tail12] = tail12-1;
	}
	head1 = tail12;

	tail3 = tail12;
}

template<typename T, template<typename T> class D>
TripleList<T, D>::~TripleList() {
	delete[] next1;
}

template<typename T, template<typename T> class D>
inline void TripleList<T, D>::removeFree(long j) {
	next1[prev1[j]] = next1[j];
	prev1[next1[j]] = prev1[j];
	if (j == freehead)
		freehead = next1[j];
}

template<typename T, template<typename T> class D>
inline long TripleList<T, D>::removeFreeHead() {
	long ret = freehead;

	next1[prev1[freehead]] = next1[freehead];
	prev1[next1[freehead]] = prev1[freehead];
	freehead = next1[freehead];

	return ret;
}

template<typename T, template<typename T> class D>
inline void TripleList<T, D>::addFreeHead(long i) {
	next1[i] = freehead;
	prev1[freehead] = i;
	prev1[i] = tail12;
	freehead = i;
}

template<typename T, template<typename T> class D>
inline long TripleList<T, D>::insertLevel2(long j, long i) {
	long h;
	if (next2[j] < tail12) {
		j = next2[j];
		while (next2[j] < tail12 && key3[vert[j]] < key3[vert[i]])
			j = next2[j];
		if (next2[j] < tail12 && key3[vert[j]] == key3[vert[i]])
			h = j;
		else if (next2[j] == tail12) {
			h = j;
			j = tail12;
		}
		else
			h = prev2[j];
	}
	else {
		h = j;
		j = tail12;
	}

	if (h == j) {
		next3[i] = vert[j];
		prev3[j] = i;
		node2[i] = j;
		prev3[i] = tail3;
		vert[j] = i;
	}
	else {
		long l = freehead;
		prev2[next2[freehead]] = prev2[freehead];
		freehead = next2[freehead];

		next2[h] = l;
		prev2[l] = h;
		next2[l] = j;
		prev2[j] = l;

		vert[l] = i;
		next3[i] = tail3;
		node2[i] = l;
		prev3[i] = tail3;
		j = l;
	}

	return j;
}

template<typename T, template<typename T> class D>
inline long TripleList<T, D>::insertLevel1Before(long k) {
	long j = removeFreeHead();
	next1[prev1[k]] = j;
	prev1[k] = j;
	if (head1 == tail12)
		head1 = j;

	return j;
}

template<typename T, template<typename T> class D>
inline long TripleList<T, D>::removeLevel1(long j) {
	next1[prev1[j]] = next1[j];
	prev1[next1[j]] = prev1[j];
	if (j == head1)
		head1 = next1[j];
	addFreeHead(j);

	return j;
}

template<typename T, template<typename T> class D>
inline long TripleList<T, D>::removeLevel2(long j) {
	if (next2[j] == tail12)
		return tail12;

	j = next2[j];
	long i = vert[j];
	vert[j] = next3[i];

	if (vert[j] == tail3) {
		next2[prev2[j]] = next2[j];
		prev2[next2[j]] = prev2[j];
		addFreeHead(j);

		if (next2[node1[j]] == tail12)
			j = removeLevel1(node1[j]);
	}

	return i;
}

template<typename T, template<typename T> class D>
inline long TripleList<T, D>::removeLevel3(long i) {
	long j = node2[i];
	long ret = node1[j];
	if (vert[j] == i) {
		if (vert[j] == tail3) {
			next2[prev2[j]] = next2[j];
			prev2[next2[j]] = prev2[j];
			addFreeHead(j);
		}
		else {
			vert[j] = next3[i];
			prev3[vert[j]] = tail3;
		}
	}
	else {
		next3[prev3[i]] = next3[i];
		prev3[next3[i]] = prev3[i];
	}

	return ret;
}

static int graphs_compar(const void* a, const void* b, void * arg) {
	return ((long *) arg)[*(long*)a] - ((long *) arg)[*(long*)b];
}

template<typename T, template<typename T> class D>
inline long ColoringHeuristic<T, D>::bucketSort(long * color, long cmin, long cmax, long * R, long r) {
	long c = cmax - cmin + 1;
	long b[c];
	long next[r+1];
	long v[r];
	long i;
	for (i = 0; i < c; i++)
		b[i] = r;
	for (i = 0; i < r; i++) {
		next[i] = b[color[i]];
		b[color[i]-cmin] = i;
		v[i] = R[i];
	}
	long j = 0,k;
	for (i = 0; j < r; i++) {
		k = b[i];
		while (k < r) {
			R[j++] = v[k];
			k = next[k];
		}
	}
}

template<typename T, template<typename T> class D>
class Greedy : public ColoringHeuristic<T, D> {
	friend class ColoringHeuristic;

	Greedy(long n, Factory<D<T> > * f) : CoveringAlgorithm<T, D>(n, f) {};
	~Greedy();

public:

	long run(Graph<D<T> > * g, long * res, long * R, long r);
};

template<typename T, template<typename T> class D>
ColoringHeuristic<T, D> * ColoringHeuristic<T, D>::greedy(long n, Factory<D> * factory) {
	return new Greedy<T, D>(n, factory);
}

template<typename T, template<typename T> class D>
long Greedy<T, D>::run(Graph<D<T> > * g, long * res, long * R, long r) {
	long i, v, co;

	//induced subgraph
	bitR->allToCompl();
	for(long j = 0 ; j < r; j++)
		bitR->toSet(R[j]);

	//initialize the colors, each vetor represents the set of colors avaible for
	//each vertex.
	for(long j = 0 ; j < n; j++)
		colors[j]->allToCompl();

	//greedy coloring
	while (i < r) {
		//in order visit each vertex of R
		v = R[i++];
		co = colors[v]->firstInCompl();
		//determine its color
		for(long j = 0; j < r; j++)
			if (R[j] == v)
				res[j] = co;

		//remove the color to the possible colors of its neigbours
		DirectlyAddressedSet::Iterator* it = g->neig(v)->newIterator(bitR, andFilter);
		long iti = it->next();
		long iu = g->neig(v)->universeSize();
		while(iti < iu) {
			colors[iti]->toSet(co);
			iti = it->next();
		}
	}

	return 0;
}

/*!
 * \class Dsatur
 *
 * Sorts input array of vertices according to saturation degree.
 *
 */
template<typename T, template<typename T> class D>
class Dsatur : public ColoringHeuristic<T, D>, private TripleList<T, D> {
	friend class ColoringHeuristic;

	Dsatur(long n, Factory<D<T> > * f) : CoveringAlgorithm<T, D>(n, f), TripleList<T, D>(n, f) {};
	~Dsatur();

public:

	long run(Graph<D<T> > * g, long * res, long * R, long r);
};

template<typename T, template<typename T> class D>
ColoringHeuristic<T, D> * ColoringHeuristic<T, D>::dsatur(long n, Factory<D> * factory) {
	return new Dsatur<T, D>(n, factory);
}

template<typename T, template<typename T> class D>
long Dsatur<T, D>::run(Graph<D<T> > * g, long * res, long * R, long r) {
	long v, u;

	//creates a bit map to the elements in R to be used as a filter
	bitR->allToCompl();
	long j;
	for (j = 0 ; j < r; j++) {
		bitR->toSet(R[j]);
		invR[R[j]] = j;
	}

	//initializes the vector of degrees
	long maxdeg = 0, mindeg = n;
	for (j = 0; j < r; j++) {
		degR->allToCompl();
		BitMapDirectlyAddressedSets<T>::intersectOf(degR, bitR, g->neig(R[j]));
		key3[j] = degR->cardOfSet();
		if (key3[j] > maxdeg)
			maxdeg = key3[j];
		if (key3[j] < mindeg)
			mindeg = key3[j];
		res[j] = j;
	}

	// sorts vertices according to degrees
	ColoringHeuristic<T, D>::bucketSort(key3, mindeg, maxdeg, res, r);
	head1 = removeFreeHead();
	next2[head1] = tail;
	key1[head1] = 0;

	v = insertDeg(head1, R[res[0]]);
	node1[v] = head1;
	for (j = 1; j < r; j++) {
		v = insertDeg(prev2[v], R[res[j]]);
		node1[v] = head1;
	}

	//initialize the vector of colors
	for (j = 0 ; j < r; j++)
		colors[j]->allToCompl();

	//coloring
	long i = r, color;
	while (--i >= 0) {
		//color the first in the heap
		v = removeDeg(head1);
		color = colors[v]->firstInCompl();
		res[invR[v]] = color;

		//update the saturation degree of its neighbors
		DirectlyAddressedSet::Iterator* it = g->neig(v)->newIterator(bitR, andFilter);
		u = it->next();
		long iu = g->neig(v)->universeSize();
		while (u < iu) {
			long w;
			if (!colors[u]->isInSet(color)) {
				colors[u]->toSet(color);
				v = removeLevel3(u);
				w = next1[v];
				long s = key1[v];
				if (next2[v] == tail12)
					removeLevel1(v);
				if (w == tail12 || key1[w] > s+1) {
					w = insertLevel1Before(w);
					next2[w] = tail12;
					key1[w] = s+1;
				}
			}
			else {
				w = prev2[node2[u]];
				removeLevel3(u);
			}
			key3[u]--;
			insertLevel2(w, u);
			u = it->next();
		}
		delete it;
	}

	return 0;
}

/*!
 * \class MCR
 *
 * Sorts input array of vertices according to saturation degree.
 *
 */
template<typename T, template<typename T> class D>
class MCR : public ColoringHeuristic<T, D> {
	friend class ColoringHeuristic;

	MCR(long n, Factory<D<T> > * f) : CoveringAlgorithm<T, D>(n, f) {};
	~MCR();

public:

	long run(Graph<D<T> > * g, long * res, long * R, long r);
};

template<typename T, template<typename T> class D>
ColoringHeuristic<T, D> * ColoringHeuristic<T, D>::mcr(long n, Factory<D> * factory) {
	return new MCR<T, D>(n, factory);
}

template<typename T, template<typename T> class D>
long MCR<T, D>::run(Graph<D<T> > * g, long * res, long * R, long r) {
	long v, u;

	//creates a bit map to the elements in R to be used as a filter
	bitR->allToCompl();
	long j;
	for (j = 0 ; j < r; j++) {
		bitR->toSet(R[j]);
		invR[R[j]] = j;
	}

	//initializes the vector of degrees
	for (j = 0; j < r; j++) {
		degR->allToCompl();
		BitMapDirectlyAddressedSets<T>::intersectOf(degR, bitR, g->neig(R[j]));
		removeFree(j);
		key1[j] = degR->cardOfSet();
	}

	//initializes the vector of degrees
	long maxdeg = 0, mindeg = n*n;
	for (j = 0; j < r; j++) {
		degR->allToCompl();
		DirectlyAddressedSet::Iterator* it = bitR->newIterator(g->neig(R[j]), andFilter);
		key3[j] = 0;
		u = it->next();
		long iu = g->neig(v)->universeSize();
		while (u < iu) {
			key3[j] += key1[u];
			u = it->next();
		}
		if (key3[j] > maxdeg)
			maxdeg = key3[j];
		if (key3[j] < mindeg)
			mindeg = key3[j];
		res[j] = j;
	}

	// sorts vertices according to degrees
	ColoringHeuristic<T, D>::bucketSort(satdeg, mindeg, maxdeg, res, r);
	head1 = removeFree(res[0]);
	next2[head1] = tail12;
	v = insertLevel2(head1, R[res[0]]);
	node1[v] = head1;
	for (j = 1; j < r; j++) {
		if ()
		v = insertDeg(prevdeg[v], R[res[j]]);
		satnode[v] = sathead;
	}

	//coloring
	long i = r;
	while (--i >= 0) {
		//color the first in the heap
		v = removeDeg(sathead);
		res[invR[v]] = i;

		//update the saturation degree of its neighbors
		DirectlyAddressedSet::Iterator* it = g->neig(v)->newIterator(bitR, andFilter);
		u = it->next();
		long iu = g->neig(v)->universeSize();
		while (u < iu) {
			long w;
			v = removeVert(u);
			w = nextsat[v];
			long s = satdeg[v];
			if (nextdeg[v] == tail)
				removeSat(v);
			if (w == tail || satdeg[w] < s-1) {
				w = insertSatBefore(w);
				nextdeg[w] = tail;
				satdeg[w] = s-1;
			}
			else {
				w = prevdeg[degnode[u]];
				removeVert(u);
			}
			degree[u]--;
			insertDeg(w, u);
			bitR->toCompl(u);
			u = it->next();
		}
		delete it;
	}

	return 0;
}

template<class T>
long ColoringHeuristic::mcrCoverSort<T>(Graph * g, long** cover, Factory<T> * factory, long* R, long r) {
	long n = g->nverts();
	long degin[n];
	long degree[n];
	long i = r;
	long v;

	//maximum minimum degree
	long maxmindeg = 0;

	//to use the same function of comparation
	graphs_deg = degree;
	graphs_sat = degin;

	memset(degin, 0, n*sizeof(long));
	memset(degree, 0, n*sizeof(long));

	//sets to represent R.
	T * bitR = factory->newInstance(n);
	T * degR = factory->newInstance(n);
	for(long j = 0 ; j < r; j++)bitR->toSet(R[j]);

	for(long j = 0; j < r; j++) {
		degR->allToSet();
		BitMapDirectlyAddressedSets<Z>::setIntersec(degR, bitR);
		BitMapDirectlyAddressedSets<Z>::setIntersec(degR, g->neig(R[j]));
		degree[R[j]] = degR->cardOfSet();
	}

	long* RR = new long[r];
	memcpy(R,RR,r*sizeof(long));
	heapify(RR, r, sizeof(long), graphs_compar);
	DirectlyAddressedSet::Iterator* it;

	//ordering
	while(i > 0){
		//determine the head of the heap
		v = RR[0];
		maxmindeg = maxmindeg < degin[v]? degin[v] : maxmindeg;
		heappoll(RR, i--, sizeof(long),graphs_compar);

		//update info
		for(long j = 0; j < r; j++)if(R[j] == v) (*cover)[j] = i;
		it = g->neig(v)->newSetIterator(bitR);
		long iti = it->next();
		long iu = g->neig(v)->universeSize();
		while(iti < iu) {
			degin[iti]++;
			iti = it->next();
		}
		heapify(RR, i, sizeof(long), graphs_compar);
	}
	return maxmindeg;
}

template<class T>
long ColoringHeuristic::mcrSort<T>(Graph * g, long** cover, Factory<T> * factory, long* R, long r, long** Ro) {
	long n = g->nverts();
	long degin[n];
	long degree[n];
	long i = r;
	long v;
	long maxmindeg = 0;
	long maxdeg = 0;

	graphs_deg = degree;
	graphs_sat = degin;

	memset(degin, 0, n*sizeof(long));
	memset(degree, 0, n*sizeof(long));

	T * bitR = factory->newInstance(n);
	T * degR = factory->newInstance(n);
	for(long j = 0 ; j < r; j++)bitR->toSet(R[j]);

	//calculate the degree in the subgraph
	for(long j = 0; j < r; j++) {
		degR->allToSet();
		BitMapDirectlyAddressedSets<Z>::setIntersec(degR, bitR);
		BitMapDirectlyAddressedSets<Z>::setIntersec(degR, g->neig(R[j]));
		degree[R[j]] = degR->cardOfSet();
		maxdeg = maxdeg < degree[R[j]]? degree[R[j]] : maxdeg;
	}

	long* RR = new long[r];
	memcpy(R,RR,r*sizeof(long));
	heapify(RR, r, sizeof(long), graphs_compar);
	DirectlyAddressedSet::Iterator* it;

	//ordering
	while(i > 0){
		(*Ro)[i] = v = RR[0];
		maxmindeg = maxmindeg < degin[v]? degin[v] : maxmindeg;
		heappoll(RR, i--, sizeof(long),graphs_compar);
		(*cover)[i] = i < maxdeg ? i : maxdeg;
		it = g->neig(v)->newIterator(bitR);
		long iti = it->next();
		long iu = g->neig(v)->universeSize();
		while(iti < iu) {
			degin[iti]++;
			iti = it->next();
		}
		heapify(RR, i, sizeof(long), graphs_compar);
	}
	return maxmindeg;
}

//>>>>>>>>>>>>>>>>>Old functions inherit from the old implementation>>>>>>>>>>>>>>>>>>>>>>>>>>.

#define graphs_declareListForAssignment(r,k) 	long _i_; \
							long _next_[r+1]; \
							for (_i_ = 0; _i_ <= r; _i_++) _next_[_i_] = r; \
							long _covhead_[k+1]; \
							for (_i_ = 0; _i_ <= k; _i_++) _covhead_[_i_] = r

#define callAssignClique(g, k, cover, i, coversize)	assignCliqueByMatrix(g, k, cover, i, &coversize, _next_, _covhead_)


static inline long graphs_assignCliqueByBitParallel(long **colors, long **adj, long k, long wordsize, long adj_bitl, long i, long mincolor, long *coversize) {
      long c=0,j=0,found=0;
      unsigned long aux=0;

      for (c = mincolor; c < (*coversize); c++) {
              for(j = 0; j < adj_bitl; j++) {
                   aux= colors[c][j] & adj[i][j];    // We are not using mask to zero vertex out of scope
                   if(aux!=0) {
                      found = 0;
                      break;
                   }
                   else found = 1;
              }
	      if(found) break;
      }
      if(c <= k) {
           colors[c][i/wordsize] |= (1 << (wordsize - (i%wordsize) -1));
           if( (c+1) > (*coversize) ) (*coversize)++;
      }

      return c + 1;
}

// sort specified vertices of a specified graph in a non-increasing order of degrees, in the
// following sense:
// in a subgraph of the specified graph induced by a set of vertices V = {V_1,V_2, ... ,V_i}, for 1 ≤ i ≤ r,
// it holds that V_i always has the maximum rank and the minimum degree in V.
// exdeg is used to break ties
// 2 vectors of r longegers are defined in the C function stack
// vertices R[j] having rank[j]=min { rank[i] | i=0,1,...,r-1 } induce a regular subgraph
// mcr order is obtained by sorting R in a non-increasing order of rank
template<class T, typename Z>
bool ColoringHeuristic<T, Z>::mcrSort(Graph<T>* g, long* d, long** ng, long* rank) {
	// degrees used locally
	long r = g->nverts();
	long locald[r];
	memcpy(locald, d, r*sizeof(long));
	// used to break ties
	long exdeg[r];
	// a negative number indicates that this parameter is not known
	memset(exdeg, -1, r*sizeof(long));

	// determine min deg
	long Rmin[r];
	long nrmin = 1;
	long i = r - 1;

	long imin = 0;
	Rmin[0] = imin;
	long degrmin = locald[imin];
	long j;
	for (j = imin + 1; j < r; j++)
		if (locald[j] < degrmin) {
			nrmin = 1;
			Rmin[0] = j;
			degrmin = locald[j];
		}
		else if (locald[j] == degrmin)
			Rmin[nrmin++] = j;

	long rmin;
	memset(rank, 0, r*sizeof(long));
	while (i >= nrmin) {
		rmin = 0;
		if (nrmin > 1) {
			for (j = 1; j < nrmin; j++) {
				if (exdeg[Rmin[rmin]] < 0) {
					exdeg[Rmin[rmin]]=0;
					long l;
					for (l = 0; l < locald[Rmin[rmin]]; l++)
						if (!rank[ng[Rmin[rmin]][l]])
							exdeg[Rmin[rmin]] += locald[ng[Rmin[rmin]][l]];
				}
				if (exdeg[Rmin[j]] < 0) {
					exdeg[Rmin[j]]=0;
					long l;
					for (l = 0; l < locald[Rmin[j]]; l++)
						if (!rank[ng[Rmin[j]][l]])
							exdeg[Rmin[j]] += locald[ng[Rmin[j]][l]];
				}
				if (exdeg[Rmin[j]] < exdeg[Rmin[rmin]])
					rmin = j;
			}
		}

		rank[Rmin[rmin]] = --i;

		while (rank[imin])
			imin++;
		for (j = 0; j < locald[Rmin[rmin]]; j++)
			if (!rank[ng[Rmin[rmin]][j]]) {
				exdeg[ng[Rmin[rmin]][j]] -= locald[Rmin[rmin]];
				locald[ng[Rmin[rmin]][j]]--;
			}

		Rmin[0] = imin;
		nrmin = 1;
		degrmin = locald[imin];
		for (j = imin + 1; j < r; j++)
			if (!rank[j])
				if (locald[j] < degrmin) {
					nrmin = 1;
					Rmin[0] = j;
					degrmin = locald[j];
				}
				else if (locald[j] == degrmin)
					Rmin[nrmin++] = j;
	}

	long maxno = 0;
	--i;
	for (j = imin; j < r; j++)
		if (!rank[j]) {
			rank[j] = i;
			maxno++;
		}

	return locald[imin] == maxno - 1;
}

// sort specified vertices of a specified graph in a non-increasing order of degrees, in the
// following sense:
// in a subgraph of the specified graph induced by a set of vertices V = {V_1,V_2, ... ,V_i}, for 1 ≤ i ≤ r,
// it holds that V_i always has the maximum rank and the minimum degree in V.
// exdeg is used to break ties; maximum rank is returned
// 3 vectors of r longegers are defined in the C function stack
// if true is returned, then vertices R[j] having rank[j]=max rank induce a clique
template<class T, typename Z>
bool ColoringHeuristic<T, Z>::mcrDegreeSort(Graph<T>* g, long * maxd, long * rank) {
	// degrees used locally
	long i, j, r = g->nverts();
	long deg[r];
	memset(deg, 0, r*sizeof(long));
	// used to break ties
	long exdeg[r];
	*maxd = 0;
	for (i = 0; i < r; i++) {
		for (j = i + 1; j < r; j++)
			if (g->hasEdge(i, j)) {
				if (++deg[i] > *maxd)
					*maxd = deg[i];
				if (++deg[j] > *maxd)
					*maxd = deg[j];
			}
		exdeg[i] = -1;
	}

	// determine min deg
	long Rmin[r];
	long nrmin = 1;
	long imin = 0;
	Rmin[0] = imin;
	long degrmin = deg[imin];
	for (j = imin + 1; j < r; j++)
		if (deg[j] < degrmin) {
			nrmin = 1;
			Rmin[0] = j;
			degrmin = deg[j];
		}
		else if (deg[j] == degrmin)
			Rmin[nrmin++] = j;

	long rmin;
	memset(rank, 0, r*sizeof(long));
	i = r;
	while (i > nrmin) {
		rmin = 0;
		if (nrmin > 1) {
			if (exdeg[Rmin[rmin]] < 0) {
				exdeg[Rmin[rmin]]=0;
				long l;
				for (l = imin; l < r; l++)
					if (rank[l] == 0 && l != Rmin[rmin] && g->hasEdge(Rmin[rmin], l))
						exdeg[Rmin[rmin]] += deg[l];
			}
			for (j = 1; j < nrmin; j++) {
				if (exdeg[Rmin[j]] < 0) {
					exdeg[Rmin[j]]=0;
					long l;
					for (l = imin; l < r; l++)
						if (rank[l] == 0 && l != Rmin[j] && g->hasEdge(Rmin[j], l))
							exdeg[Rmin[j]] += deg[l];
				}
				if (exdeg[Rmin[j]] < exdeg[Rmin[rmin]])
					rmin = j;
			}
		}

		rank[Rmin[rmin]] = i--;

		while (rank[imin] != 0)
			imin++;
		for (j = imin; j < r; j++)
			if (rank[j] == 0 && g->hasEdge(Rmin[rmin], j)) {
				exdeg[j] -= deg[Rmin[rmin]];
				deg[j]--;
			}

		Rmin[0] = imin;
		nrmin = 1;
		degrmin = deg[imin];
		for (j = imin + 1; j < r; j++)
			if (rank[j] == 0)
				if (deg[j] < degrmin) {
					nrmin = 1;
					Rmin[0] = j;
					degrmin = deg[j];
				}
				else if (deg[j] == degrmin)
					Rmin[nrmin++] = j;
	}

	long maxno = 0;
//	--i;
	for (j = imin; j < r; j++)
		if (rank[j] == 0) {
			rank[j] = i;
			maxno++;
		}

	return deg[imin] == maxno - 1;
}

static inline long graphs_assignByList(long* d, long** ng, long k, long* cover, long i, long* coversize, long* clique) {
	long j;
	// current clique starts with size 1
	cover[i] = 0;

	// find clique to include vertex as the smallest
	// not appearing in the neighborhood of i
	memset(clique, 0, ((*coversize)+1)*sizeof(long));
	for (j = 0; j < d[i]; j++)
		clique[cover[ng[i][j]]] = 1;

	long c = 1;
	while (c <= *coversize && clique[c])
		c++;
	if (c <= k) {
		if (c > *coversize)
			(*coversize)++;
		cover[i] = c;
	}

	return c;
}

template<class T, typename Z>
static inline long graphs_assignCliqueByMatrix(Graph<T>* g, long k, long* cover, long i, long* coversize, long* next, long* covhead) {
	long c = 0, r = g->nverts();
	long cc = covhead[c];
	long antcc = r;
	while (!cover[i] && c < *coversize){
		// search for cover to add
		for (; cc != covhead[c+1] && g->hasEdge(i, cc); antcc = cc, cc = next[cc]);
		if (cc == covhead[c+1])
			cover[i] = c + 1;
		else {
			c++;
			for (; cc != covhead[c]; antcc = cc, cc = next[cc]);
		}
	}
	if (c == *coversize && *coversize < k) {
		(*coversize)++;
		covhead[c] = i;
		cover[i] = c + 1;
	}
	// include i in the tail of cover list
	next[antcc] = i;
	next[i] = cc;

	return c + 1;
}

#endif /* GRAPHALGORITHM_H_ */
